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We discuss the physics of the vortex state in a d-wave superconductor, using the phenomenological 
Ginzburg-Landau theory, where many novel phenomena arise from the small admixture of the s-wave 
component induced by spatial variations in the dominant d-wave. Properties of an isolated vortex 
and of the Abrikosov vortex lattice are studied by means of analytic and numerical methods. An 
isolated vortex has a considerable structure, with four "extra" nodes in the s-wave order parameter 
Oh! symmerically placed around the core and an amplitude forming a four-lobe profile decaying as 

1/r 2 at large distances. The supercurrent and magnetic field distributions are also calculated. The 
^ \ Abrikosov lattice is in general oblique with the precise shape determined by the magnetic field and s- 

■ d mixing parameter e v . The magnetic field distribution in the Abrikosov state has two nonequivalent 

O^l ' saddle points resulting in the prediction of a double peak line shape in /iSR and NMR experiments 

as a test of a d-wave symmetry. Detailed comparison is made with existing experimental data and 
new experiments are proposed to test for the predicted effects. 
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I. INTRODUCTION 



After several years of debate there is growing agreement that the symmetry of the order parameter in the high-T c 
cuprate superconductors is not a conventional isotropic s-wave, but has a more complicated structure involving nodes 
in the gap. Recent experiments sensitive to the phase of the order parameteroH provide strong evidence for the 
d x 2_ y 2 symmetry with lines of nodes along the \k x \ = \k y \ directions. Support for the d-wave symmetry also arises 
from specific heat measurementsfj and the recent observation of a nonlinear Meisner effectu Photoemission studies,u 
Josephson interferences and c-axis Josephson tunnelingu experiments have been interpreted as being inconsistent 
with a pure d-wave order parameter. However most of these inconsistencies can be reconciled by allowing for states 
of mixed symmetryH In orthorhombic materials, such as YBCO and BiSCCO, if the dominant order parameter is 
d-wave, a small s-component will be presentcl even in a strictly uniform system. In tetragonal d-wave materials, which 
will be considered in this work, this s-wave component vanishes identically in the bulk; .however it may be nucleated 



locally by perturbations which induce spatial variations of the d-wave order parameter Ji3U e.g., by external magnetic 
fields, surfaces or impurities. In the present work we consider the vortex state of a d-wave superconductor which 
results from applying a uniform magnetic field parallel to the c-axis of the superconductor. We study the properties 
of isolated vortices and of the Abrikosov vortex lattice, both of which which differ in many aspects from those found 
in conventional superconductors, owing to the induced s-wave component. These effects will play an important role 
in transport properties of high-T c materials, which in turn are crucial for all practical applications. Understanding 
the static properties of d-wave vortices is a first important step toward the description of the more complex dynamical 
effects in the presence of transport currents, surfaces, impurities, etc. 

The pr|GJplem of an isolated vortex line in a d-wave superconductor was first studied by Soininen, Kallin and 
Berlinskyliil who considered a simple microscopic lattice model for electrons with on-site repulsion and nearest neighbor 
attraction. The resulting Bogoliubov-dc Gennes (BdG) equations were solved numerically on finite clusters to obtain 
the order parameter distribution for a single vortex. It was found that a substantial s-wave component is nucleated 
near the vortex core with opposite winding of phase relative to the d-componentJ13 and a distinct four-lobe shape pf, 
the amplitude. These results were interpreted with help of the phenomenological Ginzburg-Landau-^GL) theory,E3~til 
where the non-zero s is driven by a mixed gradient coupling to the d-component. Ren, Xu and TingtJ later attempted 
a Gorkov-type derivation of the GL theory from a continuum mean field model of d-wave superconductivity and used 
the resulting free energy to discuss the qualitative properties of a single vortex. They obtained useful asymptotic: 
expressions for the behavior of the order parameters in various regions of the vortex. Wang and MacDonaldEEl 
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investigated numerically the electronic excitations inside and outside the cores of s-wave and d-wave vortices using 
the self-consistent BdG equations. They found a distinctly different behavior of the T — quasiparticle density of 
states in the .core of the d-wave vortex compared to that in the s-wave core within the same model. Very recently 
Ichioka et alEU analyzed the structure of a d-wave vortex within the quasi-classical Eilenberger formalism- Their 
results appear to agree in every aspect with the results of GL theory presented below and in an earlier letter S3 While 
the properties of an isolated vortex are now relatively well understood, those of the vortex lattice have remained 
largely unexplored. 

Much of the work discussed above is based on a particular (effective) microscopic model of superconductivity. 
However there is presently no general agreement on the fundamental mechanism of pairing in the high-T c cuprates. 
A good alternative in such a situation is to study the phenomenological GL theory, which is based only upon general 
concepts related to symmetries of the system. Application of such theory to conventional (s-wave) superconductors 
has demonstrated its ability to predict virtually all of their phenomenological properties. In Section II, we review 
the GL theory appropriate for the d-w&ve superconductor, which involves both d-wave and an induced s-wave order 
parameter generated through the mixed gradient coupling. We discuss some of the general properties of this free 
energy and derive the corresponding GL differential equations as well as an expression for the supercurrent. In doing 
this, and throughout the entire paper, we restrict ourselves to the simple case of tetragonal symmetry, described by 
the point group D4. 

Sections III and IV are devoted to the study of a single vortex and of the Abrikosov vortex lattice. Some of the 
results described here have been previously reported in a letter S3 Here we offer a more comprehensive treatment of 
the problem, and we present a number of new results. For the single vortex we first review known analytical results 
and complement these by several new observations. We then carry out numerical integration of the GL equations for 
the single vortex geometry. In the region close to the vortex core our results confirm previous work within the BdG 
framework.til In particular we find the induced s-wave order parameter which has the expected four-lobe structure 
with minima along ±x, ±y axes and maxima along the |a;| = \y\ diagonals, and the phase winding in the opposite 
sense relative to the d-wave. Farther from the core the GL theory yields new results that were inaccessible to the 
BdG treatment due to the cluster size limitations. At a distance of several coherence lengths from the core the 
winding number of the s-wave changes from —1 to +3 resulting in four "extra" nodes in the s-wave order parameter 
symmetrically placed along the ±x, ±y axes. Analysis of the asymptotic solutions shows that these nodes are 
necessarily present in the s-component, whenever pure d-wave solutions are thermodynamically stable in the bulk. 
Our numerical work supports this conclusion. Quite generally the distribution of the d-component, as well as the 
supercurrent and the magnetic field, exhibit a four-fold anisotropy, the magnitude of which is proportional to the 
relative magnitude of s. 

As was mentioned above, the problem of the vortex lattice, which forms in magnetic fields close to the upper 
critical field id C 2, has not been previously addressed for a d-wave superconductor. In of the four- fold anisotropy 
of individual vortices one may expect that the conventional triangular Abrikosov latticeliSEj will be modified, especially 
since, even in the absence of anisotropy, the difference in the free enecgy betwecruiriangular and square lattices is 
extremely small (less than 2%). Moreover, recent neutron scatteringEil and STME3 experiments reveal an oblique 
vortex lattice in YBCO in strong magnetic fields. Ip_Section IV we solve for the structure of the vortex lattice in the 
vicinity of H C 2 . We generalize the classic Abrikoso\0 treatment to the d-wave case by first minimizing the quadratic 
part of the free energy using a Gaussian variational wave function, and then forming a periodic array of vortices from 
linear combination of these functions. We include the effects of the vector potential coupling self-consistently, thus 
improving upon our original calculation!!^ which neglected these effects. The resulting vortex lattice is found to be 
oblique, with an angle between primitive vectors ranging from 60° to 90°, depending on the strength of the mixed 
gradient coupling and magnetic field and to a lesser extent on the other parameters in the L free energy. 

In Section V we summarize our results and discuss in some detail their relevance to the existing experimental data. 
We also propose experiments that might directly test some of our predictions. 

II. GINZBURG-LANDAU THEORY OF A SUPERCONDUCTOR WITH d-WAVE PAIRING 

The Ginzburg-Landau (GL) theory for a superconductor with d x 2_ y 2 symmetry has been described by JoyntS The 
free energy density is expressed in terms of two components of the order parameter, s(r) and d(r), with appropriate 
symmetries, as follows 

/ = a s \s\ 2 + a d \d\ 2 + /?iN 4 + fa\d\ 4 + P 3 \s\ 2 \d\ 2 + p 4 (s* 2 d 2 + d* 2 s 2 ) 

+ 7 s |ns| 2 + 7d |nd| 2 + 74(n y s)*(n y d)-(n ;E s)*(n ;c d) + c.c.] +h 2 /8n. (i) 
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Here II = — iV — e*A/hc, and we assume that d is a critical order parameter, i.e., we take a s = a'(T — T s ) and 
dd = a'(T — Td) with T s < Td- The use of the same temperature derivative, a', for a s and ad is justified below. This 
also allows us to set a' = 1 in the subsequent-analysis. We assume that 0i, 02, 03, 04, 7s, 7d and j v are all positive 
as it is-suggested by a simple lattice modetj and a Gorkov-type calculation within the continuum weak coupling 
theory.113 The parameters 7 are related to the effective masses in the usual way, with 7^ = Ti 2 /2m*, for i — s,d,v. 
We shall be interested in the case when pure d-wave state is stable in the bulk in the absence o£jperturbations, i.e., 
situations when \d\ > 0, s = 0. The condition for such a state to be thermodynamically stable islill 

a d <0, 2/3 2 a s + (/3 3 -2|/3 4 |)M > 0. (2) 

With a finite d-component, the second transition temperature T s will be renormalized by the fourth order invariants. 
In particular the transition to the state with finite bulk s-wave component will occur at 

rp* _ rp (fa ~ 2|Al|) . . 

s ~ 1 s O o ir, |a \\\- L <i \p) 
2/52 ~ (P3 ~ 2|/3 4 |) 

Thus, even if the bare T s is close to Td, the true transition temperature T* may be much lower. Moreover, when 
2/32 — (03 — 2 \04 1) < 0, a second transition will never occur and we may conclude that the precise value of T s is not 
very important for the physics. I— . .— . 

There are various ways of interpreting /, some of which have been discussed by JoyntL2l and by Volovik.td Here 
we provide an interpretation in terms of nearest neighbor bond fields v(r) and h(r). These fields describe the 
superconducting pairing amplitudes on the vertical (v) and horizontal (h) bonds of the square lattice representing 
the crystalline structure of the cuprate superconductor, and arise naturally in the simple meanriicld lattice models 
of superconductivity with on site repulsion and nearest neighbor attraction between electrons .Illrc3 For tetragonal 
symmetry, the free energy /& may be written in terms of these bond fields as follows: 

fb = a {\v\ 2 + \h\ 2 ) + e(vh* + hv*) + lL {\U x h\ 2 + \U y v\ 2 ) + 7t(|II,/ 1 | 2 + \U x v\ 2 ) 

+ ic [(H x h)(U x vY + (U y h)(U y vy + c.c] + /3(H 4 + M 4 ) + h 2 /8iT. (4) 

In Eq.(|J), ao = a'(T — To), and e stabilizes the relative phase of v and h. If e is positive, then a relative phase of ir 
is stabilized, and the stable state has d-wave symmetry. If e < 0, then the quadratic terms in fb are minimized when 
v and h have the same phase, giving a state with (extended) s symmetry. The first two coefficients of the gradient 
terms and 7t involve derivatives along (e.g. H y v) and transverse (e.g. il^i;) to the bond directions. In general, 
these two coefficients will be different. The fourth order terms, proportional to (3, which are included in /{, are the 
terms which would arise in the mean field theory of XY spins. In general, a mean field theory for fermions will have 
other terms. However it is instructive to consider the consequences of these simple fourth order terms, i.e., \h\ 4 + |i>| 4 . 

The orthonormal transformation, s = (h + v)/y/2, d — (h — v)j\pl, allows us to express the coefficients of Eq.(|l|) in 
terms of the coefficients in /(,. The results are: 

a s — a — e, a<j = ao + e (5) 



01 = fa = A = P, 03 = 4/3 (6) 

7s = (iL +7t)/2 + 7c, Id = {lL + 7t)/2 - 7c, lv = {lL - 7t)/2 (7) 

The statement that the same value of a' occurs in ct s and ad is equivalent to the statement that the temperature 
derivative of e is negligible in comparison to the temperature derivative of ao- If that is not the case, then this 
approximation is not valid. In what follows we shall adopt the above approximation for computational convenience, 
but we note that it is in no way essential for the conclusions of this work, and relaxing it only leads to small quantitative 
changes. The fourth order terms, \h\ 4 + \v\ 4 , generate all of the terms in Eq.Q with-comparable magnitudes; in fact 
the resulting relative magnitudes of 0i's are very close to the weak coupling values JIj The mixed gradient term, j v , 
arises from the difference in the coefficients of the longitudinal and transverse gradient terms in the bond picture. Of 
course, this difference could be zero, but that is not expected on the basis of symmetry. 

To study the implications of the above free energy (ffl) for the structure of the isolated vortex line and the vortex 
lattice the first necessary step is to write down the field equations for the order parameters. These are obtained in 
the standard way by varying the free energy (|l|) with respect to conjugate fields d* and s*. We have 

( 7d n 2 +a d )d + lv {U 2 -U 2 x )s + 20 2 \d\ 2 d + 0z\s\ 2 d + 20 4 s 2 d* = 0, (8a) 
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( 7s n 2 + a s )s + lv {n 2 y - n 2 x )d + 2/3 1 \s 2 \s + (3 3 \d\ 2 s + 2f3 4 d 2 s* = 0. (8b) 
In a similar manner, one obtains the current density in the xy plane: 

j = S?} d *^ d) + (ridrrf ] + £|[ s * (ris) + (risrs ] 

- [s*{U x d) + (U x s)*d + c.c.] + y^ [s*(Il y d) + {U y s)*d + c.c] . (9) 

In carrying through the variational procedure it is necessary to impose appropriate boundary conditions for the 
superconductor-vacuum boundary. For our two component system these turn out to be 

n • [jjid + -y v (yU y s - xll x s)] = 0, (10a) 



n- [j s Us + lv (yn y d-xU x d)} = 0, (10b) 

where n is the unit vector normal to the surface. By combining the above two equations and comparing with the 
expression for the current density @, one can easily deduce that 

n ' j|boundary — 0, 

i.e., the normal component of supercurrent vanishes, as required on the superconductor-vacuum boundary. We also 
note that for the special case of a flat boundary along say the yz plane, conditions (|To|) acquire the simple form 
Tl x s — Tl x d — 0, which is analogous to the boundary condition in the usual one component system. 

The above set of equations constitutes a complete Ginzburg-Landau theory for a superconductor with d x 2_ y 2 pairing. 
This full theory is evidently too complicated for most practical purposes, and one must resort to approximations in 
order to obtain useful results. The rest of this paper is devoted to two such approximations valid in weak and strong 
magnetic fields. 



III. NEAR Hci- ISOLATED VORTEX LINE 



When the applied magnetic field H is close to the lower critical field, H c \ , spacing between individual vortex lines 
is large and it is sufficient to consider structure of a single isolated vortex. As mentioned in the Introduction, a 
single vortex line in a d-wave superconductor exhibits rich and rather fascinating properties that have no analogue 
in conventional superconductors with a single component order parameter. In the present section we discuss these 
properties in some detail. First we review the analytical results concerning the distribution of the order parameter, 
supercurrent and magnetic field in various regions of the vortex. Second, we carry out an explicit numerical integration 
of the GL equations for the single vortex geometry to confirm and complement these analytic solutions. 



A. Analytic solutions 



As is appropriate in the case of high-T c cuprate superconductors, we shall consider strongly type-II materials, in 
which magnetic fields vary over length scale A that is much larger than the relevant coherence length £ over which 
significant variations of the order parameter can occur. In what follows we focus only on situations where magnetic 
field is parallel to the c axis of the superconductor. 

For the problem of a single vortex line it will be convenient to work in the cylindrical gauge expressed in the usual 
polar coordinates (r, tp) , 

A = <pA(r), (11) 

with 

A(r) = - f r'h{r')dr'. (12) 
r Jo 

By adopting this particular gauge we restrict ourselves to magnetic fields h(r) = zh(r) that have no angular depen- 
dence. While this is clearly not exact for the <i-wave vortex, we shall see that quite generally the part of h that is not 
rotationally invariant is small and can thus be computed as a perturbative correction to ([ill). 
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Let us first look at the behavior of the order parameter near the center of the vortex, as r — > 0. In the relevant 
region where r <C A the magnetic field can be treated as constant, h ~ /i(0) = ho and the vector potential becomes 

Mr) = \h Q r. (13) 

For the singly quantized vortex ho can be roughly estimated by requiring that the area ~ tt\ 2 contains magnetic 
flux equal to a single flux quantum $o = hc/e*. This gives ho ~ &o/ir\ 2 . The problem is now to find simultaneous 
solutions to the two GL equations (§1) for s and d, to leading order as r — > 0. Qualitatively it is clear that at the core 
(r = 0) both d and s vanish. Moving outward from the core, the amplitude of d increases and generates nonzero s via 
the mixed gradient coupling. Around r = £<j = y/ ld/\&d\ the amplitude of d starts to level off, attaining eventually its 
bulk value do = \J I &d 1/2/32, which causes |s| to reach a maximum and then decrease to as r — > oo. This qualitative 
picture suggests that in order to study the leading order behavior we may first solve Eq.(^a| for d assuming s = 0, 
and then obtain the leading behavior of s from Eq. ( |st| ) . With this assumption Eq. ( |Sa|) becomes 

(a d + ld U 2 )d + 2p 2 \d\ 2 d = 0, (14) 

which is identical to the GL equation for the conventional one component superconductor. The asymptotic solution 
to this equation near the core is well known to beO 

d(r,i P )^(d 1 r + d 3 r 3 )e i f, (15) 

where constant d 3 is given by 

ds = -§2 [l + 27r£2V*o], (16) 

and di can be obtained by full integration of Eq. (|l4|) . Note that ordinarily only the leading dependence d(r, </?) ~ d\re %v 
is quoted; however it turns out that in our case the term d^e 1 ^ is necessary to obtain a consistent expression for 
s(r, ip). In Eq.(p7j|) the factor <J> /27r£^ dividing /i is of the order of the zero temperature upper critical field H c2 (0). 
Since we are interested in the region close to H c \ , we have ho <C H C 2 and in what follows we shall consistently neglect 
terms ~ ho/ U c i compared to unity. With this simplification Eq. (|lTJ) becomes 

d 3 ~ -di/8& (17) 

The leading behavior of s(r, tp) now can be obtained from the linearized version of Eq. (|Sb|) which reads 

(a. + 7s n 2 ) s + 7l ,(n2 -n*)d = o, (is) 

by substituting for d from Eq.(|T|). Evaluating the action of the (II 2 —II 2 ) operator in polar coordinates gives 



(HJ - E£)d(r, ¥>) = - I ) dire-* 1 + d 3 r (3e^ - e^) , (19) 

which suggests that the s-component is of the form 

s{r, <p) = Sl re~ lip + s 3 r 3 e 3ltp . (20) 
Comparing the coefficients in front of different phase factors and again neglecting terms ~ ho/H C 2(0), we obtain 

*i = I ( ^2) (21) 
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and, to the same order, S3 = 0. 

In summary, the leading order behavior of the order parameter near the core is 



d{r,ip) =dire i¥ \ (22a) 



(r, <p) 



dire 



-up 



(22b) 
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The most interesting feature of this result is the opposite winding of the s-wave component relative to d. This was first 
pointed out, by Volovikli3 based on a general symmetry argument. The solution of the form ( p2| ) was also derived by 
Ren et al&B, however the explicit form of the prefactor in the s-component is a new result of this work. Knowledge of 
this prefactor will allow us to give a simple but accurate estimate of the maximum s-wave amplitude, s max = max(|s|), 
induced in the vicinity of the vortex core. In view of the fact that far outside the core s decays algebraically with 
r (see below), such an estimate is quite important for the assessment of the relative strength of the induced s-wave 
component and the various new phenomena that its presence may lead to. The estimate is based on the assumption 
that near the core d and s rise over approximately the same length scale ~ I n particular if we assume that at 
r = £d the amplitude of d is approximately half of its bulk valueE-3 do, from Eq.(22a) we have ^ d di — do/2. Assuming 
further that |s| attains its maximum also around r = £d we arrive at the following estimate, 



d 



3 7t, 
16 a4 2 d ' 



(23) 



A similar estimate was given previously by us ,113 based on a simple argument involving the competition between the 
mixed gradient term and other second order invariants in the free energy. This argument gave the correct functional 
dependence on the GL parameters, however it missed the numerical prefactor 3/16 = 0.1875, which is important when 
investigating the quantitative properties of the above solution. Comparison to the numerical results (see the following 
subsection) shows that the above estimat e (B3[ ) is correct to within about 15%, as long as s m ax < do /4. When s ma x 
becomes larger, the asymptotic solution ([22]) is no longer justified since the condition |s| <C \d\ is violated and our 
perturbative approach starts to break down. 

A noteworthy consequence of Eq.(^) is the temperature dependence of s max near T d . If we recall that close to T d 
we have d ~ y/l - T/T d and £ d ~ l/^/l - T/T d , it follows that 



(1 - T/Td 



,3/2 



(24) 



Faster decay of the s-component near Td compared to d is a direct consequence of the fact that as a non-critical order 
parameter the former is driven by the spatial variations of the latter .t3 Thus, sufficiently close to Td, the s-component 
will always be negligible compared to d, and in many aspects a d-wave superconductor will behave very much like a 
conventional single component superconductor. Eq.(|24|) also self-consistently justifies the above perturbative solution 
of the GL equations near the core which assumes |s| to be small compared to since s max /do ~ (1 — T/Td), the 
condition |s| -C \d\ will be always fulfilled sufficiently close to Td. 

The supercurrent and the local magnetic field near the vortex core can be calculated from Eq.(^) using the order 
parameters given by Eq. (|22|) . We obtain, to leading order in r, 



is = d\ 




(25) 
(26) 



Expression (|26|) for the supercurrent shows explicitly that the s-component with opposite winding of the phase relative 
to d in fact diminishes the total supercurrent, resulting in weaker shielding of the external magnetic field compared 
to the conventional superconductor. 

We next consider the region outside the core, ^ d <C r <C A. We shall assume that in this region d has already 
reached its limiting form 



d(r, (p) = doe 1 



(27) 



Because of the condition r <C A, the magnetic field can still be treated to a reasonable approximation as constant, 
and the vector potential is thus given by Eq. ([l3|) . It is, however, easy to show that coupling to the latter can be 
ignored in this region. In particular, rewriting all the relevant operators in polar cylindrical coordinates, one can 
easily show that for d(r,ip) given by (|27j ) it holds that II 2 d(r, ip) — dor~ 2 (—l + r 2 /X 2 ) 2 . Clearly, the second term 
in the brackets (which originates from the vector potential A) can be safely ignored with respect to unity, since, by 
assumption, r/\<^ 1. With some effort, one can demonstrate that the vector potential is also negligible in the terms 

(n 2 -n 2 )d(r,^). ^ ^ 

The problem of finding the asymptotic solution outside the core region reduces to solving Eq.(pq) for s(r, <p) with d 
given by ( p7| ) and A = 0, and the additional assumptions that |s| <C \d\ and |Vs| -C |Vd|. These allow one to consider 
only the linearized equation in which the relevant terms are 
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Iv {d 2 x -tf)d+ a s s + (3 3 \d\ 2 s + 2/3 4 d 2 s* = 0. (28) 

In polar coordinates, 

(d 2 x - &l) doe* = ^ (3e 3 * - e~*) d 0) (29) 
suggesting that the solution to Eq. ( p8[ ) is of the form 

8 (r,<p) = i (he'* + he 3 *) . (30) 



Substitution in Eq.(|28|) then gives 

/l ^2 7ud0 (a s +/3 3 ^)2-4M)2' ^ 
, _ _1 w 3(a fl +^ 3 d§) + 2/34d§ 

/3 ~ 2 7 " d0 (a s +/3 3 ^) 2 -4(« 2 - ^ 

Asymptotic solution of this form was first obtained by Ren et aL0 From the knowledge of the order parameters d(r, tp) 
and s(r, tp) one can compute the corresponding distributions of the supercurrent and the magnetic field. In order to 
do this consistently, one has to include corrections ~ l/r 2 to the d-component (such as were neglected in Eq.(p7|)), as 
these are needed to insure that the continuity equation V ■ j = is satisfied. The resulting formulas can be found in 
Ref. [|TJ. 

There are two important physical consequences of Eq. (|30|) . First, the slow algebraic decay of the s-component 
outside the core region means that asymptotically in the presence of a vortex, the superconductor is not in a pure 
ci-wave state, rather there is a small s-wave admixture with angle dependent relative phase. As a result, fermionic 
excitations will be gaped in this region. As demonstrated below, only at the length scale set by the penetration depth 
is the s-component cut off exponentially and a pure d-wave state is established. 

A second interesting property of t he s- wave component can be obtained-by comparing the two solutions inside and 



outside the core. Inside the core Eq.( 22b ) implies that the winding numbc the s-component is —1. The situation 
outside the core is slightly more complex, but near Td it holds that h — — 3/i [cf. Eqs . (pl|) and ( |32"| ) in the limit 
do — » 0]. Thus the phase factor e 3ltp in Eq.(|30|) will dominate the behavior of s(r,p) and the winding number far 
outside the core will be +3. For an analytic function the winding number is a conserved topological quantity which 
can be changed only by the presence of a node. This forces us to conclude that four additional positive vortices must 
exist outside the core in the s-component S3 These "extra" vortices (or nodes) are a consequence of the topological 
constraints imposed on the relative phases of s and d by the structure of the GL equations (||) . Behavior of the phases 
6 S and 8d is schematically depicted in Fig. [l], for the two asymptotic regions as given by Eqs. (p2|), ( p7j ) and (|30"|). By 
inspection of the figure one may conclude that the four extra vortices are symmetrically placed on ±sc and ±y axes, 
since the s component apparently changes sign along these directions. Finally we note that there are no extra nodes 
in the d-component and that the total magnetic flux associated with one vortex line (consisting of 1 d-wave node and 
5 s-wave nodes) is equal to one flux quantum; there is no additional flux associated with the extra s-wave nodes. 

We have argued above that the unusual nodal structure of the d-wave vortex exists at temperatures close to Td- 
It can be shown, however, that our argument has much wider validity. It is a simple matter to demonstrate that a 
complex function of the form g(f) = ae~ ltp — be 3iV with a. b > 0, will have winding number +3 for b > a (and —1 for 



b < a). Applying this criterion to s(r, tp) given by Eq.(30) and with help of relations ( |3l| ) and (|32|), one obtains the 
following inequality 

3{a s + Ml) + Widl > {a s + f3 3 d 2 ) + 6^, (33) 

as a requirement for the winding number +3 outside the core. Upon expressing dg as |ay|/2/32 and rearranging, one 
finds that this inequality coincides with the stability condition (||). It therefore follows that for all combinations of 
GL parameters consistent with stable d-wave state, the asymptotic winding number of s outside the core is +3 and 
the non-trivial nodal structure described above exists. We may conclude that the structure of the vortex core in a 
d-wave superconductor is inherently much more complicated than that of a conventional vortex. This statement is 
valid over the entire range of temperatures in which the GL theory is applicable, in magnetic fields weak enough so 
that the vortex liae can be considered isolated. Very recently, the existence of the extra nodes has been confirmed 
by Ichioka et alrJi, who investigated the distribution of order parameters near the vortex using the quasi-classical 
Eilcnberger equations. On the other hand, however, recent numerical computations within the GL theory by Xu et 
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alxl failed to find evidence for this effect. It would be most interesting to see if evidence for the non-trivial vortex 
structure can be established in an experiment. 

Finally we shall consider the region outside the core for r>A. In this region we may still assume the asymptotic 
form ( |27| ) for d{r,<p), but we can no longer treat the magnetic field as constant. Taking into account the fact that 
•C \d\ in this region, we obtain the usual London equation for the vector potential, which in the cylindrical gauge 



reads 



V 2 A 



The asymptotic solution to this equation for r A is 



A = 



2ir\ 



7T A 

2r 



1/2 



,-r/A 



(34) 



(35) 



which gives the usual exponentially decaying magnetic field far from the vortex. Using the vector potential given 
by (p5|) one can solve for s(r, <p) from Eq.(Sb). The result, to the leading order in (r/A), is 



s(r, ip) 



with 



91 



7T A 

2r 



■S3 



1/2 



-r/A 



. ' " (axe-** + 83*"*) , 



Jvd 



2A 2 



(ft - 2f3 4 )d 2 ■ 



(36) 



(37) 



Thus, as expected, the s-wave component will be exponentially small beyond distances from the core in excess of A, 
and on these large length scales the d-wave superconductor will behave as a conventional single component type-II 
material. Eqs.(^6|) and also show that to leading order, the total winding number of s(r, ip) remains undetermined 
(see the discussion of winding above). However, upon computing higher order terms in (r/A) one finds that the winding 
number in this region remains +3, so that there are no additional nodes present in the s-component. 



B. Numerical results 



The analytic results presented in the preceding subsection establish rich and complex structure of the vortex line 
in a d-wave superconductor; however owing to the rather complex structure of the underlying GL equations (|J) the 
analytic treatment is restricted to limiting cases where certain small parameters can be identified. Consequently, the 
information such a treatment provides is mainly of qualitative nature. In order to study the problem in more detail, 
we have integrated the GL equations numerically. Besides confirming the above analytic predictions, the numerical 
approach is capable of addressing the behavior of the order parameter at length scales comparable to £<j, where the 
analytic approach is difficult. In particular we will be most interested in the detailed behavior of the s-component 
near the core, focusing on its exotic nodal structure that was predicted by topological arguments. 

In order to arrive at a truly selfconsistent numerical solution, one should in principle complement the GL Eqs.(^) 
by the Maxwell equation V x h s — ^Lj an d include the induced magnetic field h s in the total vector potential A. 
However, as we are mostly interested in the region near the core (r <C A), it is permissible to neglect these screening 
effects and indeed the coupling to the vector potential altogether, provided that we impose correct boundary conditions 
for a single vortex geometry (see below). Neglecting the vector potential leads to a significant simplification of the 
problem. Physically this corresponds to the extreme type-II limit, A/£ — > oo. For a realistic system where A/£ is finite 
(but large), ignoring the vector potential coupling is equivalent to neglecting terms ~ (r/A) 2 compared to unity [see 
discussion following Eq. (p7|)] . 

With the vector potential absent from the GL equations it is convenient to introduce a set of dimensionless quantities 
such that a s is measured in units of \atd\, /3-parameters in units of 2/32,s and d in units of the bulk d-wave gap do, 
and all the lengths in units of This allows one to write the GL Eqs.^8|) in the following simple dimensionless form 

- (V 2 + 1) d + e v {d 2 x - d 2 y )s + \d\ 2 d + (3 3 \s\ 2 d + 2(3 4 s 2 d* = 0, (38a) 
- (V 2 - a.) 9 + e v {d 2 x - d 2 y )d + 2(3 1 \s 2 \s + f3 3 \d\ 2 s + 2(3 A d 2 s* = 0, (38b) 
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where e v — jv/jd and we have set 7 S = jd- On physical grounds [cf. Eq.(fj])] we do not expect 7 S and jd to differ 
dramatically; and we have verified that allowing 7 S ^ jd does not have a significant effect on the solutions. 

We have integrated Eqs.(^8|) numerically on a rectangular N x N domain for the boundary conditions appropriate 
for a single vortex: 

d\boundary d^€ ^ i ^boundary 0- (^^) 

We used-an iterative Newton's algorithm as described in Ref. [ ^8|. At each step of iteration the Conjugate Gradient 
Methodca was use to solve the resulting system of linear equations. 

Results of our numerical analysis indeed confirm all of the qualitative features found by the analytic considerations 
of the preceding subsection. Fig. ^ shows the behavior of the d- and s-wave amplitudes near the center of the vortex, 
with parameters described in the figure caption. The resulting amplitude of the s-component for this particular 
parameter configuration was s max = 0.024do, in good agreement with estimate (|2^). A domain size of N = 201 was 
used in the numerical integration, encompassing a physical region of the size L ~ 20£<2- In Fig. |^ only the central 
(121 x 121) region is displayed, where the boundary effects are expected to be strongly suppressed (the numerical 
solution was in fact well behaved all the way to the boundary of the system). As expected for this relatively weak 
admixture of the s-component, the amplitude of d has almost perfect circular symmetry. The amplitude of s is nearly 
circular in the inner core of th e vortex and it shows marked fourfold anisotropy outside the core, in accordance with 



the asymptotic solutions ( 22b ) and (|30|). Four symmetrically placed maxima along diagonals and four nodes along 
±x and ±y are visible in the contour plot. To see these more clearly we show in Fig. ^ the amplitudes of the s-wave 
component along x-axis and a x = y diagonal. A node close to x = 3£<j is clearly visible, which is nothing else than 
one of the four extra vortices. The figure also confirms the linear behavior of |s| and \d\ near the origin and the fact 
that both rise on approximately same length scale £<j. One can also see the 1/r 2 decay of |s| outside the core region, 
where |d| is constant. 

Fig. U shows the superconducting phases 9d and 8 S of the two components of the order parameter. While the 
distribution of 9d looks very much like that of conventional singly quantized vortex, the distribution of 9 S is clearly 
more complicated. In particular the opposite winding of the phase near the core and four positive vortices along the 
±x, zty axes are clearly distinguishable. Comparison to Fig. |lj shows that our numerical results are again in complete 
agreement with the analytical predictions summarized in the preceding subsection. 

The important quantity that determines the nature of excitations in the vicinity of the vortex line is the relative 
phase A8 = 6 S — 9d- We plot A8 in Fig. [|. Over much of the region the relative phase is AO = ±7r/2, resulting 
in a d ± is state that has minimum gap equal to |s|. This is a direct consequence of the fact that /3 4 > in the 
free energy (|l|). However, the phase difference cannot be equal to ±7r/2 over the entire area since this would be 
incompatible with the topological constraints that require opposite winding of the two components near the core. 
Thus narrow domain walls appear along tho-dzx, ±y axes, where AO changes rapidly. This result is in agreement with 
the microscopic treatment of Soinincn et aZilil within the Bogoliubov-de Gennes theory. However since the complexity 
of this formalism did not allow one to extend the calculations to sufficient distances from the core, the extra vortices 
were originally not found. The present GL theory, being inherently simpler allows us to study larger clusters. As one 
moves further out from the core, domain walls abruptly end at the cores of the four s-wave vortices and AO starts to 
vary more slowly while being still locked to ±7r/2 over large areas. 

Supcrcurrent j produced by the above order parameter distribution, computed numerically from Eq. ([)[), is shown 
in Fig. |^. Panel (a) shows the distribution of the magnitude | j | while panel (b) displays streamlines of the vector field 
j. Note that because of the Ampere's law V x h s = — j, the latter is equivalent to the lines of constant magnetic 
field given by the supercurrent, and thus Fig. ^(b) also gives the distribution of the spatially varying component of 
the screening field. 

In addition to the particular case described above we have numerically studied a large number of other parameter 
combinations. All show similar behavior. The feature that changes between different configurations is the relative 
magnitude of the s-component, which is, as we have explicitly verified, well described by Eq.(|2^). The larger the 
ratio s m ax /do, the more anisotropic ci-wave component becomes and along with it the distribution of supercurrent 
and induced magnetic field. As an example of such a case we show amplitudes of s and d in Fig. [?], for the particular 
parameter combination that yields s rnax ~ 0.15<io- The relevant supercurrent distributions is plotted in Fig. ^|. 

IV. NEAR H C 2- STRUCTURE OF THE VORTEX LATTICE 

In what -follows we present our treatment of the vortex lattice problem. In general we follow the path outlined by 
AbrikosoVjtj with necessary modifications that arise from the presence of two order parameters in the free energy. 
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A. Linearized GL equations and their variational solution 



In the vicinity of the upper critical field H C 2 the amplitudes of orderparameters are small, and the essential physics 
is contained in the linearized field equations that are obtained from (BH by neglecting the nonlinear terms: 



( 7d n 2 + a d )d + lv (U 2 y - U 2 x )s = 0, 



(40a) 



( 7s n 2 + a s ) s + 7 „(n2-n2)d = o. (40b) 

Formally this corresponds to the expansion to leading order in the small parameter (H C 2 — H)/H c2 . The gauge 
invariant gradient II can be separated into two pieces, 

n = ri + it = (-iV - e*A /ch) - e*A s /ch, (41) 

where H = V x Ao corresponds to the uniform applied field, and h s = V x A s is the screening field induced by the 
supercurrent j s in the sample, given by the Maxwell equation 



_ , 47T. 

V x h s = — j s 

c 



(42) 



Let us for a moment ignore complications arising from the screening effects and consider only the vector potential 
Ao. This is permissible, since as it will become clear later, corrections to Eqs.(|40|) from the screening field are of 
the same higher order in the small parameter (H C 2 — H) / H^,a.s the nonlinear terms which have been neglected in 
these equations. In the same spirit as the original Abrikosovt^l treatment, these higher order terms will be included 
variationally in a later stage of the calculation. 



It is easily seen that in the Landau gauge Ao — yHx the linearized field equations (40) are satisfied by taking 

d(r) = e lky d(x), s(r) = e iky s(x). (43) 

Thus, exactly as in the one component case first studied by Abrikosov,0 we are left with one dimensional problem 
which can be stated as follows 



P 



1 



f , -„, 2/ \2 

— + -muj c {x-x k ) 
2m 2 



P 



1 



f i * 2/ \2 



s = 0, 



(44a) 



P , 1 2/ n2 

2m 2 



P , 1 2/ N2 

— + -muj c (x - x k ) 
2m 2 



0. 



(44b) 



Here we have denoted p — —ihd/dx, x k — kl 2 , and lj c = (e*H/mc). The magnetic length I = ^/Tic/e*H will play the 
role of characteristic length for the vortex lattice. We also assume henceforth that m* d — m* = to, i.e., that jd = Js, 
and we use e v = jv/ls — m * s / m * v . Equations (Q) resemble those for the quantum mechanical harmonic oscillator 
problem with the potential centered at x — x k . In view of the fact that x k is arbitrary, it is clear that Eqs.(|44|) will 
have infinitely many degenerate solutions which can be labeled by continuous index k. This degeneracy will play a 
crucial role later when we construct the periodic space-filling solution. However, for the moment, we shall ignore 
this issue and investigate Eqs.(|44|) with x k fixed. The essential difference from the one component case is that these 
equations have no obvious exact solutions. In what follows we shall seek suitable variational solutions to Eqs.(Q). In 
order to stress the analogy with the harmonic oscillator, we may write (J44|) in the following way 



(Ho + a d )d + Vs = Ed, 



(45a) 



Vd + {Ho + a s )s = Es, (45b) 

where Tto = fiu) c (a'a + 1/2) and V = e v (hLu c /2)(a^a^ + aa) are expressed in terms of the usual raising and lowering 
operators, which can be written as a = [(x — x k )/l + l(d/dx)]/V%- By including the right hand side of Eqs.(|4§) we 
are considering a slightly more general problem: E — corresponds to the physical solution for H = H C 2(T), and 
solutions for E < will be useful later when we consider the stability of various vortex lattice structures. 
In order to motivate our variational solution to the linearized problem, let us define 



10 



H ± =H ±V, i P ± =d±s. 



(46) 



In terms of these new variables, the set of equations (|45| ) becomes 

H++T-T* -AT/2 \ ftfi + \ ( V + 

\=e[ ). (47! 

-AT/2 H-+T-T* J \ip~ ) \y~ 

where we have defined 

T* = (T d + T s )/2, AT — T d — T s , (48) 

for convenience in later calculations. A nice feature of the representation (E^) is that for the degenerate case AT = 
the equations for tp + and ip~ decouple, each becoming a simple harmonic oscillator problem. Motivated by this fact 
we shall seek the variational solution for the general case in the form of normalized ground state wave-functions of 
the harmonic oscillator, 



^ {x ) = ^ e -°l^-^l*\ (49) 
The variational parameters a + and <r_ will be determined by minimization of the eigenvalueE] 

(E) = (T — T*) + ^+H + tp+) + \{<p-H-<p-) - ~AT<^-), (50) 

where angular brackets stand for spatial averages. All the integrals are easily evaluated and if one defines er+ = 
crcos??, <7_ = crsini?, the resulting expression for (E) can be explicitly minimized with respect to a 2 . The minimum 
occurs for a 1 — tani? + 1/ tant?, and is 



(E) _T-T* If huj. 



AT AT 4 V AT 



(l + e^)tan?? + (l-e^) ' 



tan?? 



1 / 2 tain? 

2VTTta^' (51) 



The last equation must be minimized numerically with respect to tan??. It is also clear from this equation that two 
parameters, e v and 

A = hu c /AT, (52) 

determine the nature of the variational solution. In the two limiting cases the exact minimum can be easily found. In 
the low field limit, A <C 1, we have cr+ w er_ w 1, while in the high field limit, A ^> 1, we have a± w [(l±e 1 ,)/(l=Fe„)] 1 / 4 . 
It follows that at least intermediate values of A are required for appreciable effects from s-d mixing to occur. Otherwise 
(p+ ~ (p— anc [ according to Eq.(^6|) the s-component effectively vanishes. 

Solutions to Eq.(|l|) with (E) = give the dependence of the upper critical field H c i on the temperature. Whenever a 
finite admixture of the s-component is present, a characteristic upward curvature is found near the critical temperature 
in H C 2(T). Such curvature has been observed |6sperimentally in both LSCO and YBCO compounds!^ and has 
been interpreted as a consequence of s-d mixing. OEJ For given parameters T d and T s and several values of e v such 
dependence is shown in Fig.^, as obtained by numerical minimization of Eq. (pl|) . 

B. Vortex lattice solution 

To construct a periodic vortex lattice consider a linear superpositions of the basis functions (|4^) of the form 

$ ± (r)=^ Crl e m 'V±(,), (53) 

n 

where c„ are complex constants. In order to impose periodicity in y-direction we have constrained the values of k to 
integer multiples 

k n = nq, n = 0,±l,±2,... (54) 
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of the parameter q which will be determined later from the requirement of minimum free energy. The space filling 
solutions of GL equations can be written as, cf. Eq. (|46|) , 

d(r) = [¥+(r)+tf_(r)]/2, 

s(r) = [* + (r)-*_(r)]/2. (55) 

These solutions will also be periodic in x pipyided that the constants c„ satisfy the condition c n+ M = c n for some 
integer N. As was first noted by Abrikosov,tj the analysis of the vortex lattice for general N is extremely difficult. 
It was however conjectured that the absolute minimum of the free energy takes place for N < 2, in which case the 
analysis is simplified. In what follows we shall restrict ourselves to the case of TV = 2 for the two component system. 
Taking N = 2 we have c 2n = Co and c 2n +\ = c\. This, along with Eqs.(|54|) and (|55"|), implies periodicity of s and d in 
x and y with periods 

L x = 2l 2 q, L y = 2ixjq. (56) 

Each rectangular L x x L y unit cell then contains an amount of flux 

HL x L y = 2(hc/e*) = 2$ (57) 

where $o stands for the flux quantum. Thus, by construction, each rectangular unit cell contains exactly two singly 
quantized vortices, independent of the value of parameter q. The resulting vortex lattice may be thought of as centered 
rectangular with two quanta per unit cell or, equivalently, as an oblique cell with lattice vectors of equal length and 
one flux quantum. While the restriction to centered rectangular lattices is made primarilv-foii the computational 
convenience, it is also compatible with recent experiments on YBCO which show evidence ! 2 HP 2 ! for oblique vortex 
lattices with nearly equal lattice constants in high fields. 

The parameter q controls the shape of the unit cell. It is customary to define the ratio 

R = L x /L y = {l 2 /n)q 2 1 (58) 

and it follows that R = 1 corresponds to the square, R = \fi corresponds to the triangular, and the intermediate 
values 1 < R < V3 imply the oblique vortex lattice. 

The solution that we have constructed for the GL equations (JsJ) has three free parameters, Co, c\ and R. These 
parameters determine the structure of the vortex lattice near H C 2- Within the linearized approximation to the GL free 
energy these solutions are degenerate in energy. It is the fourth order terms that lift this degeneracy and determine 
the equilibrium lattice structure. In order to find this minimum one must take into account the fourth order terms in 
the free energy (|l|) as well as the effects of screening which were so far ignored. 

The complete average free energy density (Q) can be written as 

</) = (fa) + (fa) + {h 2 )/&n, (59) 

where fa and fa stand for quadratic and quartic invariants respectively, and h = H + h s is the local magnetic field. 
Let us now consider the effect of screening by looking at the gradient terms in (fa) with the complete LI as given by 
Eq. ([ill) ■ A typical term will be of the form 

(|rid| 2 ) = (\u d + u s d\ 2 ) ~ (|iM 2 ) + (n s • [d*u d + c.c.]), (60) 

where in the last equality terms of the order of |n s (i| 2 have been neglected. This is consistent with the general idea of 
GL theory of keeping only terms up to fourth order in the order parameters. Being proportional to the supercurrent, 
LI S already contains terms quadratic in the order parameters. If we expand all the remaining gradient terms in the 
similar way, systematically neglecting terms containing order parameters to powers higher than four, we can write 
the result as 

(/ 2 ) = (/ 2 (0) ) + 4(n s -j s ). (6i) 

e 

Here is the part of fa containing only the IIo piece of the gauge invariant gradient, i.e., the quadratic part in 
the absence of screening, and similarly j s is assumed to be given by Eq.(|^) with n = n . If we take into account the 

property of the variational solution (f^) = E(\s\ 2 + \d\ 2 } that follows from Eqs.(45) and use the definition of n s we 
can write 
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(/> = E(\s\ 2 + M| 2 ) - (l/c)(A s • j s > + (f 4 ) + ((H + h s ) 2 )/8^, (62) 

The second term on the RHS can be simplified by expressing j s through the Maxwell equation (fl2|). Integrating by 
parts and neglecting the surface term one obtains (l/c)(A s ■ j s ) = (h 2 )/An. 

Similarly the last term on the RHS can be rewritten recalling the definition B = H + (h s ) of the magnetic induction 
as (/i 2 )/8tt - H 2 /8n + B • H/4tt. 

The manipulations performed above are useful since in fixed applied magnetic field the proper thermodynamic 
potential to minimize is the mean Gibbs free energy density related to / by (g) = (/) — B • H/Air. For this quantity 
we finally arrive at an expression 

(g) = E(\s\ 2 + \d\ 2 ) + </ 4 ) - (h 2 s )/87r H 2 /8tt. (63) 

Before we proceed with minimization of the Gibbs potential let us notice that the simple thermodynamic relation 
d(g)/dH = —B/4tt can be used to extract the average screening field in the superconductor 

(h.) = B-H = -(^j{\a\ 3 + \d\ 3 ). (64) 

A similar relation between the average induced field and the order parameter for the conventional s-wave 
superconductor£L3 is known as the "first Abrikosov identity," but the corresponding determination of the spatial 
distribution of h s (r) is more complicated (see below). It is easy to verify that in the limit e v — > (i.e. in the limit of 
pure (i-wave) Eq.(|64|) assumes the precise form of this identity, including all the relevant prefactors that follow upon 
expressing dE/dH from Eq.(|l]). In the vicinity of H c2 it holds that E ~ (8E/dH)(H - H c2 ) and it follows that to 

the leading order we can write (f^) — (h s )(H C 2 — H)/A-k. This allows us to express the Gibbs free energy in the 
form where the leading dependence on the magnetic field H is manifestly displayed: 

(g) - (g)n = ]-(h C 2 - H)(h s ) + (U) - ^(hl), (65) 

with (g) n — —H 2 /8n being the normal state contribution to the Gibbs free energy. 

Consider now a simple scaling transformation (s, d) — > (As, Ac?) where A is a real number. It is clear that under such 
a transformation (h s ) — > X 2 (h s ), while (fi) — > A 4 (7.4) and (h 2 ) — > X 4 (h 2 ). Consequently, the Gibbs free energy ( |65| ) 
will have a well defined minimum for the particular value of A. We use this property to determine the normalization 
of the order parameters s and d. Carrying out the minimization we obtain 

.9 - (.9 n = -g 3-777 77727-, 66) 

8vr 87r(/ 4 ) - (hi) 

an expression which is clearly independent of the particular normalization of s and d. If we further define the Abrikosov 
ratio Pa and the Ginzburg-Landau parameter k by 

a ^s) 2 A (h) la >7\ 

^ = jh^' =^WY (67) 

we can write the resulting Gibbs free energy for the Abrikosov vortex lattice in the familiar fornB 

Several remarks are in order. The Abrikosov ratio (3a defined by Eq.(|67]) is independent of the coefficients fa in the 
quartic part of the free energy and depends only on the shape of the unit cell in the vortex lattice. To the extent 
that k is independent of the specific lattice shape, the minimum Gibbs free energy corresponds to the minimum 
of /3a, which generalizes the familiar Abrikosov result, apart from writing it in terms of magnetic field instead of 
the absolute squared order parameter. As will be shown below by numerical calculation, it is indeed true that the 
parameter k defined by E q.(|67| ) depends only very weakly on the vortex lattice shape, and thus the factor (2k 2 — 1) 
in the denominator of Eq. (|68[) serves simply as the criterion for type-II behavior, which occurs only for k > l/y/2. It 
is in this sense that one can think of k as a generalization of the conventional Ginzburg-Landau parameter; we note 
that k defined by Eq.(^) cannot be simply related to the usual ratio of penetration depth A to coherence length £. 
This difficulty is related to the fact that in the presence of two order parameters s and d we have, strictly speaking, 
two distinct coherence lengths, £ s and Most observable phenomena will only reveal a single "effective" coherence 
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length given by a certain combination of £ s and but this will presumably depend on the type of probe used in 
the experiment. By contrast, there will be only single penetration depth A, as this quantity is related to the decay of 
magnetic field inside the superconductor. Alternatively, A may be viewed as a measure of the bulk superfiuid density, 
which is in the present case associated with the d-wave component only, since the s-wave vanishes in the bulk. Thus 
it may be suggested that k — \/£a, where £a is the effective coherence length relevant to the Abrikosov lattice, 
determined by the usual criterion of overlapping vortex cores at H = H C 2- 



C. Magnetic field distribution 



The ultimate goal of this section is to determine the actual structure of the vortex lattice by minimizing the Gibbs 
free energy given by Eq. (|68|) . To obtain the parameters (3 a and k that enter this expression it will be necessary 
to evaluate the spatial averages (fa) and (h 2 s ) (note that quantity (h s ) has been already calculated in Eq.(|64|) by 
thermodynamic argument). The former of the two averages can be computed in a fairly straightforward manner since 
fa is directly related to the vortex lattice solutions \l/±(r), which are simple linear superpositions of the Gaussian 
wavefunctions <p k given by Eq.(^£). The situation with the other average, (h^), is more complicated as one has to 
first invert the Maxwell equation (42[) in order to express the local screening field h s (r) in terms of the supercurrent 
j s . Both of these quantities are themselves of interest, as they can be in principle measured by various experimental 
probes (see Sec. VI for the discussion). 

With this in mind let us calculate the spatial distribution of the screening field. If we express h s in terms of the 
vector potential A s , the Maxwell equation (E^) can be written as 

Air 

V 2 A S = j„ (69) 

c 

where we have taken advantage of the fact that the Landau gauge satisfies V • A s = 0. The easiest way to invert the 
equation ([39]) is to exploit the periodicity of the vortex lattice solution and work in Fourier space. If we write 

j.(r) = Xy k - r j s (k), A s (r) = £e ik ' r A s (k), (70) 

k k 

where the summation goes over the reciprocal lattice vectors k = (k Xl k y ) = (2irki/ L x , 2-nkijLy) and ikxik^) = 
0, ±1, ±2, . . ., Eq.@ implies that 

A*(k) = ^#, MO. (71) 
Thus, one obtains for the Fourier components of the screening field 

h s (k)^^, MO. (72) 

In order to evaluate this expression it is helpful to write the supercurrent (|^) using wavefunctions \&± instead of s 
and d: 

j,(r) = 8 e — [x(l - ae„)**n x * Q + y(l + ae v )^* a U y ^ a + c.c] . (73) 

c*=± 

In order to model the layered structure of cuprate superconductors we have introduced the usual geometrical factor 
5 = (layer thickness/layer spacing). The case 5 — 1 corresponds to the cubic lattice, while (5^0 represents the limit 
of a single isolated layer. In this notation, the Fourier components of the supercurrent are 

* Pi 

j.(k) = 6— J2 - ae v )(^* a U x ^ a ) k + y(l + ae v )(** a ILy* a ) k + (c.c.)_ k ] . (74) 

Q = ± 

Here we have introduced a shorthand notation 

1 ("La fLy 



L X Ly JQ 



(...) k = / dx dy...e~ lkr , (75) 
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which will prove very convenient in the subsequent calculations. With some effort, the following useful relations can 
be derived: 



2 \<7 a 



ik y a a (|* tt | 2 )k, 



(*;n,* Q ) k = 



Za a \<T a 



(76) 



In the real space, ^(r) is a linear combination of Gaussians, and thus the Fourier components (| x I'a| 2 )k are easily 
evaluated. One obtains 



«k 



T~ ex P i -T[(^/ a «) + (*w°"a) ] r ' 



(77) 



where 



a k = e 



c 4 2 + (-i) fel ci4 2+1 ] 



(78) 



are real constants, independent of particular lattice shape. Substituting relations ( |76| ) in the expression for the 
supercurrent (173) one obtains 



j.(k) = ^ [£(1 - ae u )(-fc, CT 2 ) + y(l + ae,)(fc x /4)] <|* Q | 2 ) k . 



(79) 



a = ± 



This expression is particularly useful for numerical evaluation of the real space supercurrent distribution, since in 
view of the Gaussian dependence of (|4 , a | 2 )k on k [cf. Eq.(|77j)1 it is clear that the corresponding Fourier series will 
converge very rapidly. 

Finally, we are in the position to give the local screening field. Substitution of the above equation ((79[) into the 
Maxwell equation ( |72| ) yields all the Fourier components of the field with k ^ 0. The k = component is just the 
real space average of the screening field (h s ) given by Eq.(Q). Combining these results we obtain, after some algebra, 
the real space field distribution of the form 



h s (r) = — zwd- 



Zo £ e lk - r (|* + | 2 + |*_| 2 ) k + ^5> 

k k#0 



k r ( K K 

\ k 2 + k% 



(l*+| 2 -|*-| 2 > 



(80) 



where we have defined numerical factors 



z =[(a 2 _+a%)+e v (a 2 _-al)]/2, 
Zl =[(a 2 _-a 2 + )+e v (a 2 _+a 2 + )}/2. 

We notice that the first Fourier sum in the brackets of Eq.Q is equal to |* + (r)| 2 + |*_(r)| 2 = 2(|s(r)| 2 + \d(r)\ 2 ). 
Thus in the limit of pure <i-wave state where e v — > and a± — * 1 the correspondence with the Abrikosov result for 
a conventional superconductor becomes transparent. In this limit we have zq — > 1, z\ — > and |s(r)| — > 0, and the 
spatially varying form of the Abrikosov first identity is recovered, with d(r) playing the role of the conventional order 
parameter ^E'(r). The second sum clearly has a nonlocal dependence on the order parameters and can be written as 
/ d 2 r'g(r — r')[s(r')e?*(r') + s* (r')d(r')]. Such a term has no counterpart in the conventional theory, and arises only 
as a result of mixing between s- and d-components of the order parameter. The nonlocality of this term is a direct 
consequence of the symmetry of the problem, since by itself the term (sd* + s*d) is not invariant under D^, it can 
enter only in combination with other terms of proper symmetry. 



D. Structure of the vortex lattice 



As mentioned above, in order to determine the shape of the vortex lattice, one needs to evaluate the averages of 
fourth order terms {fi) and (h 2 ). Now that the distribution of the magnetic field h(r) has been derived, evaluation 
of these averages is a straightforward, albeit quite a lengthy procedure. The technical details of this calculation are 
worked out in the Appendix, and here we only summarize the results and discuss some of the physical implications. 
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Equations (A6) and (All) of the Appendix give the expressions for the fourth order averages (fi) and {h s ) in terms 
of rapidly converging sums that are suitable for numerical evaluation. Making use of these, the Abrikosov ratio and 
Ginzburg-Landau parameter can be expressed in the following simple way 

P A = §H' [(2 ° + z ^k)^++(k) + (z - z 1 ry k )fi__(k)] 2 , (82) 

k 

and 

K 2 = 4E k Eq[Mi»L(k) +M2»aa(k)» aa (k) +/i 3 » act (k)n aa (k) +/x 4 nL(k)] (83) 
irS 2 (e*h/mc) 2 J2kl( z o + ^i??k)^++(k) + (z - ^i?7k)^ — (k)] 2 

where the prime on the sums means that only terms with k\ and k 2 both even or odd are included, fl al 3(k) are 



Gaussian functions given explicitly by Eq.(A5), r?k is a simple function defined by Eq.(A9, and a = —a. Note that 
the above expression for Pa is independent of parameters Pi that enter the quartic part of the free energy /4, and 
other quadratic parameters enter only via the variational parameters <r±. 

The shape of the vortex lattice unit cell is determined by the ratio R — L x /L y . The value of R that corresponds to 
the thermodynamically stable configuration, R m im is obtained by requiring that the Gibbs free energy is minimum. 
Equation (^) shows that, at given external magnetic field H, the Gibbs free energy (g) is entirely determined by the 
two parameters given above, Pa and n. Numerical evaluation of these parameters confirms that n is only very weakly 
dependent on the particular lattice shape, as it is illustrated by Fig. [lO[ The dependence of the Gibbs free energy ( |6"g| ) 
on R is almost entirely contained in the Abrikosov ratio [3a and thus, in most of the parameter space, the minimum 
of (g) coincides to a good accuracy with the minimum of (3a- For example in the particular case displayed in the Fig. 
|l0| the minimum of [3a differs by less than 2% from the minimum the full free energy. 

Fig. ^ also shows a typical dependence of (3a on R for different values of mixed gradient coupling e v , as obtained 



by numerical evaluation of Eq.(82). When e v — the superconductor is in a pure d-wave state with no s-wave 
component present. Within the phcnomcnological GL theory, this situation is identical to the case of a conventional 
superconductor studied by Abrikosov. Thus, the state with minimum free energy has R m in = Vo which corresponds 
to thp-,usual triangular vortex lattice. In this limit we obtain the correct value of (3a — 1.1596 as quoted by Kleiner 
et alEQ However, as soon as non-zero coupling e v is introduced, the situation changes and the minimum of (3a shifts 
to the values R m in < V3, signalling that an oblique vortex lattice is favored. The minimum R m in varies continuously 
with e v and at certain value of £„, which depends on other parameters in GL free energy, R m in reaches the value of 
1, corresponding to the square lattice. Further increase of e v then has no effect on the shape of the lattice, which 
remains square. 

One may conclude that in a d-wave superconductor, in the regime close to the upper critical field H c2 , a general 
oblique vortex lattice is thermodynamically stable, unless the material is in one of thc-.limiting regimes in which 
the mixed gradient coupling e v is very small or very large. NumericaO and analyticalEj calculations based on the 
simple mean field model with proper symmetries, find evidence for a mixed gradient term of about the same order of 
magnitude as the conventional gradient terms. This would seem to argue against the two limiting cases mentioned 
above. 

An example of the oblique vortex lattice is displayed in Fig. [0], where we show the d and s components of the order 
parameter as obtained by numerical evaluation of Eqs.(p5|), for given set of GL parameters. An interesting conclusion 
can be drawn by comparing the two components of the order parameter: it is evident that the non-trivial nodal 
structure of the s-wave component, such as was described in Sec. Ill for an isolated vortex, persists in this high field 
regime. Indeed, zeros of s are present in the regions where \d\ > 0. This quite remarkable result appears to suggest 
that the "extra" vortices in the s-component are present over the entire portion of the phase diagram representing 
the mixed state of a d-wave superconductor. 

Many experimental probes are sensitive to the spatial variations of the magnetic field rather than to the order 
parameter itself. The spatially varying component of the magnetic field, h s (r), as evaluated from Eq.(^) is shown 
in Fig. [l^. Notice that as a consequence of the Maxwell equation V x h s = (47r/c)j s , it follows that the contours 
of constant magnetic field coincide with the supercurrent streamlines. Comparison to the order parameter plot in 
Fig. |ll| confirms that magnetic field and supercurrent distributions have the same symmetry as the vortex lattice. 
A non-trivial nodal structure of the s-wave has an effect on the field distribution, which develops two nonequivalent 
saddle points, marked SI and S2 in Fig. |lj. In principle, it might be possible to determine such structure by /iSR or 
NMR experiments. Fig. [lj| displays the /iSR line shapes that result from the magnetic field distribution as discussed 
above. The quantity shown is 

P{h) = -L- [ 6[h - h(r)]d 2 r, (84) 

J-> X ljy J 
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for the case of triangular, oblique and square flux lattices. In the triangular and square lattices symmetry requires 
only one type of saddle point, resulting in the conventional single peak structure. In the oblique lattice, which is 
characteristic of a d-wave superconductor, the two non-equivalent saddle points give rise to two distinct Van Hove 
type singularities. Appearance of two distinct peaks in /iSR or NMR spectra, would provide evidence for d-,W|ave 
behavior, since the explanations of oblique vortex lattice that invoke anisotropy within single component modelE3 do 
not lead to this effect. 

The last subject that we want to address here concerns the spatial orientation of the vortex lattice with respect to 
the crystalline axes of the superconductor. From Fig. [H], it can be seen that the principal axes of the vortex lattice 
are not aligned with any of the high symmetry directions of the underlying crystal. Instead, it is the (110) direction of 
the vortex lattice that coincides with the (100) or (010) directions of the crystal. It turns out that the construction of 
the vortex lattice as presented above forces this particular orientation and does not allow for identical configurations 
that are rotated by some angle a. In the traditional one component case, this is not a concern since the free energy 
has full rotational invariance. In the present case, however, we must take a closer look at these rotated configurations 
as we have terms in the free energy that break rotational invariance. It is conceivable that such rotated configurations 
might in fact be lower in free energy than the ones we have considered so far. In what follows we show by an explicit 
calculation that this is not the case, and that we have in fact found the solution that corresponds to the absolute 
minimum of / as given by Eq.(|l|). 

Consider a simple rotation of the coordinate system in the xy plane by an angle a 

x — x cos a — y sin a, 

y = x sin a + y cos a. (85) 

Under such transformation all the second order terms are invariant except for the mixed gradient term which trans- 
forms as follows 

ds_dd*_ _ ds_dd*_ + _ ( co „2 Q sip 2 x f ds dd* ds dd* 

dy dy dx dx \dy' dy' dx' dx' 

f ds dd* ds dd* \ 
+2 S1 nacosa^— — + — — j+c.c. (86) 

One can now derive and analyze the linearized field equations using the rotated coordinates (x',y') in exactly the 
same way as before, and the angle a becomes just another variational parameter with respect to which the free energy 
is minimized. It turns out that it is possible to write down the linearized equations for s and d that are identical to 
Eqs.(E^) with V changed to V — e v (hui c /2)(e 2 ' ia aJaJ + e~ 2la aa). In such case, one expects there will be a constant 
phase difference 2a between s and d components, and the appropriate variational solution is of the form 

d(x) = e- ia [<p+(x) + V -(x)}, s(x) = e ia [<p+(x)- i p-(x)}, (87) 

where (p + and <p~ are the normalized lowest eigenfunctions of harmonic oscillator as defined by Eq.(^). The energy 
eigenvalue is easily evaluated, and we obtain a generalization of Eq, 

(E) _T-T* 1 (HluA 

' hair \ 1 — tan 2 d 



1 / 2tan$ 



2 V 1 + tan 2 



(1 + e v cos 2a) tanz? + (1 - e v cos 2a) 
1 + e„ sin 2 2a 



AT J 1 + tan 2 



It is a simple matter to minimize (E) with respect to a, and one finds that minima can occur only for a — 0, ±7r/2, ±7r. 
Thus, we are led to the conclusion that within our variational solution the most stable vortex lattice is the one aligned 
with the underlying crystal as described above (cf. Fig. |TT| ) . Among the fourth order terms in the free energy only 
f3 4 (s* 2 d 2 + s 2 d* 2 ) depends on a. This dependence is particularly simple; upon rotation the constant (3^ changes to 
/34COs4a. Clearly, this term only has minima for trivial values of a = 0, ±7r/2, . . ., so the above conclusion should 
hold even when the fourth order terms are included. In order to verify that this conclusion is not altered by some 
complicated interplay between angular dependencies of /a and fi, we have carried out the numerical minimization of 
the free energy of the rotated vortex lattice, along the lines indicated for the case a = 0. We find that, for all the 
regions of the parameter space that were investigated, the absolute minimum of the free energy occurs for a = 0. As 
a consistency check we have also verified that identical minima are found for a = ±7r/2,±7r, which corresponds to 
the discrete rotations under the D4 group. 

The above conclusion concerning the orientation of the vortex lattice may be understood by analyzing the mixed 
gradient term in the free energy density ([!]). Its structure forces the vortex lattice to align in such a way that the 
greatest gradient of order parameters is along one axis, while the smallest possible gradient is along other axis. An 
arrangement of vortices such as the one shown in Figure [ll] definitely satisfies this requirement. 
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V. SUMMARY AND DISCUSSION 



The main goal of this work was to present a detailed study of the vortex state in a d-wave superconductor, focusing 
on the properties arising from s-d mixing that have no analogue in conventional superconductors. Analysis of the 
vortex state is done in two regimes; in the vicinity of H c \ where the properties of isolated vortex lines can be studied, 
and near H C 2 where the collective properties of vortices forming a dense lattice are important. 

For the single vortex line the topological structure of the induced s-wave order parameter is highly non-trivial, 
consisting of one counter-rotating unit vortex, centered at the core, surrounded by four additional positive vortices 
located symmetrically at a distance of several coherence lengths from the core. A new result of this work is the 
realization that the above structure will occur for all parameter configurations that give rise to stable d-wave in the; 
bulk (provided one is well into the type-II regime), and not only in the vicinity of T c as was originally suggested.ES 
This conclusion is confirmed by an explicit integration of the GL equations over the wide range of parameters, 
and also by the calculations of Ichioka et aZO, who find analogous behavior using the quasi-classical Eilenberger 
equations. The question arises as to whether this non-trivial topological structure of a single vortex could be probed 
experimentally. There are clearly many complicating factors which are likely to render this task very difficult. The 
main challenge arises from the fact that one expects the induced s-component to be small, on the order of few percent 
of d. Such a small admixture of s might be hard to detect directly, and the corresponding distortion of the d-wave 
amplitude, supercurrent and magnetic field distributions will also be small. It might in principle be possible to probe 
the s-component by scanning Josephson tunneling from an s-wave tip, which by symmetry would not couple to the 
dominant d-wave. With sufficient resolution such an experiment could detect strong anisotropy in the s-component 
and possibly also the extra nodes. The internal structure of a vortex will also have an effect on the transport properties; 
e.g. it is conceivable that it may lead to changes in the Magnus force acting on a vortex in a current field. These 
issues clearly require further investigation. 

Finally we note that although a finite induced s-component will restore the gap along the \k x \ = \k y \ directions in 
the vicinity of the core, this will not invalidate the prediction of VolovikS regarding the ~ y/~H contribution to the 
density of states (DOS) on the Fermi surface, which was recently confirmed by specific heat measurement by Moler 
et ala Volovik's prediction is based on the observation (originally used by Yip and SaulscJ to predict the nonlinear 
Meissner effect) that the quasiparticle excitation spectrum is shifted by the superfluid velocity field around the vortex 
core, with the dominant contribution coming from quasiparticles far from the core in position space and close to the 
nodes in momentum space. Since the amplitude of the s-component far from the core vanishes as 1/r 2 the reduction 
of the DOS will be always negligible beyond a certain distance from the core compared to the energy shift due to 
superfluid velocity which decays only as 1/r. Thus at relatively small fields compared to H C 2, such as were used in the 
specific heat measurements, there will be no correction to the Volovik's result from the induced s-wave. At stronger 
fields, when the vortices are closely spaced, corrections may appear; however, in such a case one expects Volovik's 
derivation to break down since the concept of an isolated vortex with a well defined asymptotic flow field loses its 
meaning in the dense Abrikosov lattice. 

The vortex lattice near H C 2 is in general oblique for a d-wave superconductor. The precise shape determined by 
an angle <fi between primitive vectors depends in a complicated way on the parameters in the GL free energy; most 
strongly on the mixed gradient coupling e v and on magnetic field via the parameter A = Huj c / AT, which also determine 
the relative magnitude of s. Quite generally, when e„ = 0, the s-component vanishes and the lattice is triangular. 
By increasing e v and A the lattice is continuously deformed, becomins-oblique and eventually square. Observation 
of an oblique flux lattice with </> ~ 73° was reported by Keimer et alBH using small angle neutron scattering from 
YBCO in magnetic field*i0.5 T< H < 5 T. This would be in agreement with our result, although as was pointed out 
by Walker and TimuskJ23 the observed distortion may also be accounted for by the intrinsic a-b plane anisotropy of 
the orthorhombic YBCO crystal. More recently an oblique vortex lattice with <fi ~ 77° was found in YBCO using 
STM by Maggio-Aprile et <d£3 This technique also revealed elongated vortex cores with the ratio of principal axes 
about 1.5. If, as noted by authors, this elongation was due to the a-b anisotropy in coherence lengths, within a 
simple London model of s-wave superconductivity this would lead to the flux lattice with an angle inconsistent with 
the actual observed value of 77°. Thus it would appear that the a-b anisotropy alone cannot explain the observed 
distortion in the vortex lattice and additional effects, such as the internal symmetry of the order parameter, must be 
invoked in order to account for the experimental data. In this respect it would be most interesting to see if an oblique 
lattice can also be established experimentally in truly tetragonal superconductors. Alternatively it would be desirable 
to study the analogous GL theory for the Z?2 orthorhombic symmetry; unfortunately such a theory is complicated 
and contains even more phenomenological parameters so that a quantitative comparison with experiment would be 
difficulties An alternative way of distinguishing between the effects of a-b anisotropy and d-wave symmetry is to study 
the magnetic field distributions in the vortex lattice. The present theory predicts a double-peak structure in uSR or 
NMR lincshapes whenever the flux lattice is oblique, while interpretations based on simple scaling argumentsES lead 
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to conventional single-peak lineshapes. 
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APPENDIX: EVALUATION OF QUARTIC AVERAGES </ 4 ) AND (h 2 s ) 



We first evaluate the contribution of (fi). For the purposes of calculation it is convenient to express (fi) in terms 
of the functions x &±, 



(h) = Ml*+| 4 ) + m 2 (I*+I 2 I*-I 2 } + M3(l*+l 2 *+*-> + M^*- 2 } + 

where the constants /i are given as follows 
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(A2) 



The easiest way to evaluate the spatial averages is to express them as Fourier series. For example, one can write the 
typical member as follows, 



(A3) 



where we have used only the basic properties of the Fourier series. The utility of this formulation lies in the fact that 
components of the form ($* can be expressed in terms of simple Gaussians, and consequently the summations 
indicated in Eq. (|A3|) converge very rapidly. In particular it is useful to define 



(*a*/3>k = a k rj a/3 (k), 



(A4) 



where the coefficients a k are given by Eq. (|78|) . The factors f2 a ^(k) contain all the dependence on the lattice structure 
and can be evaluated by explicit integration; we have 



n aa {k) = -J-exp|--|[(/c a; /cr a ) 2 + (k y (7 a ) 2 ]j , 

H Q/3 (k) = £ V /=P^ exp {-?-£-r[kl + hi + aj)k x k v }} , a^0. 

In terms of these functions, (/4) can be written in a compact form, 

(/ 4 ) = ^ a k a -k X! lMi^ 2 a( k ) + M2^aa(k)ri aa (k) + /i 3 £l QQ (k)ft aa (k) + /i 4 0^(k)], 

k a=± 



(A5) 



(A6) 



which is suitable for numerical evaluation. In deriving Eq. ( |A6| ) we have used the symmetry Q| g(k) = Q a p(— k) which 
is apparent from Eqs.(A5), and we use the notation a — — a. 

Let us now turn to calculation of (h 2 s ) . A similar approach as above will work here if we write 



= J>.(k)ft.(-k). 



(A7) 
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The Fourier components h s (k) can be deduced from Eq. (| 



ft.(k) = -nS^ [z (|*+| 2 + |*-| 2 )k + ^k(|*+| 2 - |*-| 2 )k] , (A8) 
where we have introduced the quantity 

k 2 -k 2 

»fr=< kl+k * T (A9) 

0, if k = 

which allows all of the Fourier components of h s (k) to be expressed by a single equation ( |A8| ) . A compact expression 
for the average squared field can be written in terms of the functions Q ai g(k) and has the following form 



(h 2 s )=ir 2 5 2 l Va k a_ k [(z +^i?7k)n ++ (k) + (zo-2i?7k)0__ (k)] 2 . (A10) 

V mc I * — ' 

N ' k 

Let us finally write, for the purpose of completeness, the expression for the average field: 

(h s ) = h s (k = 0) = -Tc5z ao (—) [n++(0) + n__(0)] = -2nS^- (—) (|c | 2 + |d| 2 ). (All) 

\ mc J L x \ mc J 

We have now expressed all the averages that are needed to minimize the Gibbs fre e ener gy (pq ) . In principle, both 
Pa and k can now be evaluated numecically using the rapidly converging sums (A3) and ( |A10| )7 However, following 
the work of Kleiner, Roth and Autler,c3 it is practical to carry out the minimization with respect to the constants cq 
and ci analytically. This reduces the relevant parameter space to a single parameter, the geometric ratio R — L x /L y . 
To this end we have purposely singled out the dependence on these parameters in the expressions for (/ 4 ) and (h 2 ). 
In particular, the entire dependence on cq and c\ is contained in the constants ak, and both (f4) and {h 2 } are of the 
form 

^oko-k/^.Aa) = £ {e^ fc ^[ Co 4 2 + (-i) fe i Cl 4 2+1 ]} 2 /(fc 1 ,fc 2 ), (A12) 

k ki k2 

where f(ki, fe) is independent on cq and c%. Our goal here is to factor out the entire dependence of this expression on 
the constants cq and C\. This is done by considering separately the cases when k\, hi are even and odd, and recalling 
that by assumption C2k = Co and c 2 k+i = C\. We obtain an expression of the following form 

(M 4 + |ci| 4 ) t/( 2fc i> 2fc 2) + /(2fci + 1, 2*a)] 

+ (cocf + c* 2 c 2 ) J2 [/(2fci, 2*2 + 1) - /(2*i + 1, 2fc 2 + 1)] (A13) 

+ 2|c | 2 | Cl | 2 [/(2*i, 2fe) - /(2*i + 1, 2fc 2 ) + /(2fe l5 2* 2 + 1) + /(2*! + 1, 2fc 2 + 1]. 
fcifc 2 

Clearly, the entire denominator 8'7r(fi) — (h 2 ) of the free energy ( |66| ) can be written in the a bove form, where the 
dependence on constants Cq and C\ is explicitly shown. Combining this with the expression ( All ) for the average 
induced field (h s ) it is easy to see that the total free energy (^) can be written schematically as 

(M 2 + H 2 ) 2 (A14) 

(\c \^ + \c 1 \^Go(R)+2\c \ 2 \c 1 \ 2 G 1 (R) + 2Re(c 2 c* 1 2 )G 2 (R)' [ 1 

where Gi(R) are complicated functions of R and other GL parameters, but are independent of Cq and c\. The 
above expression can be easily minimized with respe ct to c and c\ ; one obtains that a condition for the minimum is 
ci = ±icq. Clearly, the value of the expression (A14) only depends on the ratio Ci/cq, so we can arbitrarily choose 

Co = 1, c\ = i. (A15) 

With this choice, we have an identity 
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a k a_ k = [(-l) fcl + (-l) fc ^, (A16) 



which will simplify evaluation of the sums in {fn) and (h 2 s ). Also, the expression for the average induced field (All) 
simplifies: 

(h s ) = -4kS^- ( —) . (A17) 



L r V mc 
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FIG. 1. Schematic diagram of phases 6d and 9 a of the two components of the order parameter in the two asymptotic regions 
close to and far from the vortex core. Note that this diagram is more complete than the similar one published in Ref. [ hll in 
that it includes the region outside the core. The present diagram also differs from that in Ref. [ [l5| which shows (we believe 
incorrectly) the s-component with opposite overall sign outside the core. 

FIG. 2. Contour plot of the amplitudes of the (a) d-wave and (b) s-wave components of the order parameter. The GL 
parameters used for the plot are: 7 S = 7^ = 7„, a 3 = 10|ad|, Pi = Pa = 0, and Pi = 0.5/32- The lightest regions of the diagram 
correspond to the largest amplitudes. The scale is in units of the d-wave coherence length £<i. 

FIG. 3. Amplitude of the d-wave and the s-wave component along the z-axis (solid line) and along the diagonal x = y 
(dotted line) normalized to the bulk value do. The parameters used are the same as in Fig. ^| The d-component is almost 
completely isotropic for this case so that the two cuts are indistinguishable. 

FIG. 4. The angle of arrow with respect to the horizontal x-sods represents the phase of the (a) d-wave and (b) s-wave 
components of the order parameter. The parameters used are the same as in Fig. ^. 

FIG. 5. Relative phase A8 of the two components of the order parameter, for the same parameters as in Fig. |^. 

FIG. 6. Contour plot of (a) supercurrent amplitude (b) supercurrent streamlines (which coincide with the lines of constant 
magnetic field), for the same parameters as in Fig. ^ 

FIG. 7. Contour plot of the amplitudes of (a) d-wave and (b) s-wave components of the order parameter for a different set 
of GL parameters: 7 3 = j d = 7«, ct s = 1.4ja d j, Pi = P2 = Ps = Pi- 

FIG. 8. Contour plot of (a) supercurrent amplitude (b) supercurrent streamlines (which coincide with the lines of constant 
magnetic field), for the same parameters as in Fig. |^. 

FIG. 9. Dependence of upper critical field H C 2(T) on temperature for various values of parameter e v and T s = 0.5T d , 
T = 0.75T d . 

FIG. 10. Abrikosov ratio Pa as a function of the lattice geometry factor R — L x /L y for different values of e v . Note that the 
minimum of Pa is moving from R = \/3 to R — lase„ increases. This implies a continuous deformation of the initially triangular 
vortex lattice into an oblique and square lattice. The parameters used are T s = 0.5T d , T — 0.75T d , P\ = P2 = Pz = Pi = 1, and 
B — 0.8H C 2- The inset shows the /^-dependence of the squared Ginzburg-Landau ratio k 2 on approximately thesame scale. 

FIG. 11. Contour plot of the amplitudes of (a) d-component and (b) s-component of the order parameter in the vortex 
lattice. The same parameters are used as in Fig. holwith e„ = 0.45 resulting in an oblique vortex lattice with R m in = 1-29 and 
the angle between primitive vectors <f> = 76° . The oblique unit cell containing one flux quantum is marked by a solid line. 

FIG. 12. Distribution of the magnetic field h s in the vortex lattice. Letters M, m, SI and 52 denote the maximum, 
minimum and two saddle points respectively. GL parameters used for the plot are same as in Fig. [ll| 

FIG. 13. Typical /xSR lineshapes as obtained from the magnetic field distribution in the vortex lattice. The curves shown 
are for triangular (e„ = 0.0), oblique (e„ = 0.4) and square (e v = 0.6) flux lattices. Magnetic field on the horizontal axis is in 
the units of ho = \{h 3 )\ = 4TrS(zo/L x )(e*h/mc) and the curves are offset verticaly for clarity. 
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